load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"

begin
    a = addfile("/public/users/wrf/ncepdata/Reanalysis/2011/pres.sfc.2010.nc", "r")
;    printVarSummary(a)

    ; read time
    time = a->time
    ntime = dimsizes(time)
    ; convert UTC time to CST time
    time = time + 8
    timeStr = ut_string(time, "%Y-%N-%D_%H")
    timeStr_output = ut_string(time, "%Y-%N-%D_%H:%M:%S")

    ; read pressure
    pres = a->pres
    lon = a->lon
    lat = a->lat
;    printVarSummary(pres)
    presNew = (/pres*pres@scale_factor+pres@add_offset/) * 0.01 ; Unit: hPa
    copy_VarAtts(pres, presNew)
    copy_VarCoords(pres, presNew)
;    presNew!0 = "time"
;    presNew!1 = "latitude"
;    presNew!2 = "longitude"
;    presNew&latitude = fspan(-90., 90., 73)
;    presNew&latitude@units = "degrees_north"
;    presNew&longitude = fspan(0., 357.5, 144)
;    presNew&longitude@units = "degrees_east"
    delete(presNew@scale_factor)
    delete(presNew@add_offset)
;    printVarSummary(presNew)

    ; plotting
;    wks_type = "x11"
    wks_type = "png"
    wks_type@wkHeight = 1200
    wks_type@wkWidth = 1200

    res                     = True
    res@cnFillOn            = True
    res@gsnSpreadColors     = True
    res@gsnAddCyclic        = False
    res@mpDataBaseVersion       = "MediumRes"
    res@mpDataSetName           = "Earth..4"
    res@mpGeophysicalLineColor      = "Black"
    res@mpNationalLineColor         = "Black"
    res@mpUSStateLineColor          = "Black"
    res@mpGridLineColor             = "Black"
    res@mpLimbLineColor             = "Black"
    res@mpPerimLineColor            = "Black"
    res@mpGeophysicalLineThicknessF = 2.0
    res@mpGridLineThicknessF        = 1.0
    res@mpLimbLineThicknessF        = 1.0
    res@mpNationalLineThicknessF    = 2.0
    res@mpUSStateLineThicknessF     = 1.0
    res@mpOutlineBoundarySets       = (/"National","Continental", "Geophysical"/)   ; turn on country boundaries
    res@mpOutlineSpecifiers = "China:states"

    ; crop China Region
    minLat = 10
    maxLat = 50
    minLon = 70
    maxLon = 130
    pres_China = presNew(:, {minLat:maxLat}, {minLon:maxLon})
    res@mpMinLatF           = minLat
    res@mpMaxLatF           = maxLat
    res@mpMinLonF           = minLon
    res@mpMaxLonF           = maxLon
    ;    res@cnLevelSelectionMode = "ManualLevels"
    ;    res@cnLevelSpacingF     = 40.
    ;    res@cnMaxLevelValF      = 1100.
    ;    res@cnMinLevelValF      = 600.
    ;    res@lbAutoLabelStride   = True


    ; for Gansu North Grid
    pres_Gansu_North = presNew(:, {40.7140}, {96.2402})
    printVarSummary(pres_Gansu_North)

    ; for Gansu South Grid
    pres_Gansu_South = presNew(:, {34.3525}, {104.8754})
    printVarSummary(pres_Gansu_South)

    ; output Gansu Pressure to TXT file
    line_array = new(ntime + 1, "string")
    l_cnt = 0

    line_array(l_cnt) = "Time, Gansu North(40, 95), Gansu South(35, 105)"
    l_cnt = l_cnt + 1
    do it = 0, ntime-1, 1
        line_array(l_cnt) = timeStr_output(it) + \
            sprintf("%10.4f", pres_Gansu_North(it)) + \
            sprintf("%10.4f", pres_Gansu_South(it))
        l_cnt = l_cnt + 1
    end do
    pres_txt_filename = "~/Data/NCEP/pres_2010.txt"
    asciiwrite(pres_txt_filename, line_array)

    do it = 0, ntime-1, 1
        output_path = "~/Data/NCEP/figures/sfcPressure/"
        system("mkdir -p " + output_path)
        output_filename = output_path + "NCEP_sfcpres_" + timeStr(it)
        wks = gsn_open_wks(wks_type, output_filename)
        gsn_define_colormap(wks, "BlAqGrYeOrRe")

        pres_China_it = pres_China(it, :, :)
        pres_China_it@long_name = "Surface Pressure at " + timeStr(it)
;        printVarSummary(pres_China_it)
;        printVarSummary(pres_China)
        plot = gsn_csm_contour_map(wks, pres_China_it, res)
;        print(max(presNew(it, :, :)))
;        print(min(presNew(it, :, :)))
    end do

end


